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We model the behavior of N classical impurity fields immersed in a larger Bose-Einstein condensate by 
N+l coupled nonlinear Schrodinger equations in 1, 2, and 3 space dimensions. We discuss the stability of the 
uniform miscible system and show the importance of surface tension for self localization of the impurity fields. 
We derive analytically the attractive tail of impurity-impurity interaction due to mediation by the underlying 
condensate. Assuming all impurity fields interact with the same strength, we explore numerically the resulting 
phase diagram, which contains four phases: /) all fields are miscible; //) the impurity fields are miscible with 
each other but phase separate from the condensate as a single bubble; IIP) the localized impurity fields stay 
miscible with the condensate, but not with each other; and IV) the impurity fields phase separate from the 
condensate and each other, forming a crystalline structure within a bubble. Thus, we show that a crystal can 
be constructed solely from superfluid components. Finally, we argue that the crystalline phases maintain their 
superfluid behavior, i.e. they possess a nonclassical rotational inertia, which, combined with lattice order, is a 
characteristic of supersolidity. 

I. INTRODUCTION 

Bose-Einstein condensates (BECs) in trapped atomic gases iH 0], being highly manipulable and well approximated by a 
classical nonlinear mean field theory, have proven to be ideal systems in which to realize many manifestations of nonlinear 
physics, such as bright and dark solitons, quantized vortices and their lattice formation, modulation instabilities, and so forth 
(see iH] for a review). In this paper, we use these same qualities of dilute BECs to explore other aspects of nonlinear physics, 
namely the self localization, phase separation, and crystallization of impurity fields embedded within a larger condensate. 

Specifically, we investigate the formation of nontrivial impurity structures, such as a crystal composed solely of superfluid 
components, by modeling the impurity fields and condensate by coupled nonlinear Schrodinger equations (NLSEs). We are 
therefore dealing with a type of multicomponent condensate mixture — a system that has long generated much interest (see for 
example chapters 15 and 16 of JH] and references therein) — composed of one large component and many smaller components 
that we will refer to as impurity fields. The NLSEs, acting as a nontrivial yet often tractable model of an idealized superfluid 
system, have a long and distinguished history of providing insights into the fundamental nature of superfluidity. Indeed, when 
dilute Bose-Einstein condensates in trapped atomic gases were discovered, the NLSEs were often used successfully to describe 
quantitively these new superfluid systems. 

The analysis in this paper provides the first step toward a theoretical understanding of realistic experiments where distin- 
guishable coherent impurity fields in trapped gases might be achieved by employing different atomic levels, isotopes, species of 
atoms, or some combination thereof. This analysis builds upon the rich history of investigating impurities in Bose fluids, which 
have been widelyused not only for probing the properties of Bose fluids, but also for creating new phases of matter (see for 
example S B B H B E lll^ Some of the results presented in this paper were first reported in lll2|] . 

In the following section, we introduce the model considered in this paper — specifically, + 1 coupled NLSEs — and 
point out the relevant conserved laws. The criteria for the system to collapse and the instability criteria of uniform miscible 
state are put forward in section III. With these established, we proceed to study the existence of a single localized structure as 
a solution of (1+1) coupled NLSEs in section HV] Here, we outline variational arguments for the existence of this solution for 
various dimensions, examine the importance of surface tension, and compare it with direct numerical simulations. In section IVl 
we develop a perturbative expansion to derive the condensate-mediated interaction among impurity fields. Section IVll presents 
the different phases that arise in this system and we show that crystallization of the impurity fields is possible in two different 
regimes. We then focus on the case of an impurity crystal forming within the bubble immersed in a condensate in section Ivnl 
Finally, in section IVIIII we show that the condensate field displays nonclassical rotational inertia in the crystal regime and, as 
such, bears some of the primary features of supersolidity, namely lattice ordering and nonclassical rotational inertia. 
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II. MODEL 

In this paper, we consider a large Bose-Einstein condensate denoted by ijj coupled to N small distinguishable impurity fields 
(treated as classical fields) denoted by x/e- We make the assumption that all impurity fields interact with the same coupling 
constants, which allows for an uncluttered description of the system's nontrivial properties. Our system of TV + 1 coupled 
classical fields is therefore governed by the following Hamiltonian: 



H 




N N N N \ 

k=l /c=l /c=l j<k J 



(1) 



7o is the self interaction of the impurity fields, 7 is the interaction between impurity fields, and A is the coupling of the condensate 
to the impurity fields. Unless otherwise stated, we shall restrict ourselves to positive parameters. In ultracold dilute atomic 
condensates these coupling constants are directly proportional to the atomic scattering length with a proportionality constant 
of Anh^ /m. As mentioned above, in principle, these scattering lengths can be tuned in current ultracold trapped atomic gas 
experiments. 

Using this Hamiltonian, the dynamics of the system are governed by + 1 coupled NLSEs (in the condensate context the 
equations are often referred to as Gross-Pitaevskii equations 111 [2]): 



1 ^ 



\Xk\ 



N 



(2) 
(3) 



where A stands for the Laplace operator in D spatial dimensions. In this system there is conservation of the mass of the large 
condensate field (particle number) A^*^^^ = / \ilj\'^d^x and that of each impurity field Uk = J \xk\'^d^x, and we assume that 
Uk <^ N^^\ The total energy eq. ([B and the total linear momentum 



r ^ I r 

Im / rVV^d^x + Im^- xl^Xkd'' 



of the system are also conserved. Note that the momentum of each individual field is not conserved. 



III. COLLAPSE OF THE SYSTEM AND INSTABILITIES OF THE UNIFORM STATE 



In this paper, we are interested in the nontrivial structures that impurity fields can generate in a condensate. To this end, first 
note that there are two regimes to avoid. One is the finite-time collapse of the system, and the other is the uninteresting uniform 
miscible state. To find out when these undesirable states occur, consider the (A^ + 1) x (A" + 1) interaction matrix 



M 



/ 1 A A 

A 70 7 

A 7 70 



A A \ 

7 7 

7 



(4) 



70 7 
7 70 / 



Mp where the vector of the densities is given by p 



A : 
VA 7 . 

derived from the quadratic form of the potential energy, i.e. 

Ixip7 IX2p, •••)• One can show that if M is negative semidefinite the system will experience a finite-time collapse ill 
Furthermore, because the density is non-negative, one can extend the criteria in [14] and show that the system would also 
experience a finite-time collapse if M is conegative or positive subdefinite defined as a ^ • ^MiuXiXk < for all > (for 
a discussion on the properties of these matrices see lllSll ). We avoid situations that are proven to exhibit finite-time collapse by 
our assumption that the coupling constants are positive. 

Next, by ignoring kinetic energy terms (which means ignoring surface tension and effects due to the system being of finite 
size), one can prove that the uniform state is energetically stable if and only if M is positive semidefinite, i.e. if all eigenvalues 
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of A4 are non-negative It is when a modulational instability arises that the nontrivial structures emerge. We will therefore 
focus on situations where at least one eigenvalue of M is negative. 

The + 1 eigenvalues of M. are 70 — 7, which is N — 1 degenerate, and 

I ((1 + (A/" - 1)7 + 70) ± + - 1)7 + 70)^ + 4A/'A2) . Therefore, when 7 > 70 or A > V(^^7+7o ^ 

system becomes energetically unstable. This highlights the difference between the two types of instabilities in this system, 
which is important for the analysis in this paper. In the first regime, i.e. 7 > 70, each of the impurity fields localizes with respect 
to the other impurity fields regardless of whether or not they are phase separated from the condensate. In the second regime, i.e. 

f A/"— 1) 7+7o 

A > ^ , the large background condensate phase separates from the impurity fields. For a discussion of the dynamical 

instabilities of the miscible state, see Appendix lAl 



IV. SELF LOCALIZATION OF A SINGLE IMPURITY FIELD 



In this section, we will discuss the system of a single small impurity field embedded in and interacting with a large condensate. 
We will show, by a variational argument, that there is a critical coupling parameter between the condensate and impurity beyond 
which the impurity field self localizes. By self localization we mean the ground state of the system will contain an impurity field 
with a finite extent (or with an extent smaller than the size of the system) due to the interaction between the impurity field and 
the surrounding condensate. We will show that surface tension allows one to interpolate between the previous known results of 
self localization of a single impurity atom [11] and the bulk phase separation condition [16]. 

We estimate the ground-state wave functions via a variational argument for the system with a single impurity field in 1, 2, and 
3 dimensions. We impose that the trial wave functions possess the following characteristics: that x(^) be localized such that 

D_ 

x(r) ^ as r ^ 00, and its particle number be fixed such that = Cd J \x\'^^^~^dr^ where Cd = is the surface of a 

unitary sphere in D spatial dimensions. The condensate wave function should be depleted where the impurity field is positioned, 
and ilj{r) ipo = constant as r ^ 00. (Note ipo is real because it is a ground state.) The appropriate energy to be minimized 
is 

l|w./.i2 I 1/L/,|2 I/, \2\2 I \/L/,|2 i/, |2\i i2 , ^ |V7.j2 , 70, 



E = CdJ (^2l^V'r+2(IV'l -|V'o|^)^+A(|^|^-|Vonixr + ^|Vxl^ + f Ixrir^-^rfr. (5) 
We will consider the following real, normalized trial functions: 




X{r) = ^^^^ /(q^O (6) 

^{r) = ^o{l-ax{rf). (7) 

with variational parameters a and a > 0, both determined once the energy estimates are minimized (note that, in a real and 
finite system, the limit a ^ cannot exist); and /(r) is localized, i.e. /(r) as r oo. The normalization constant is 
JVd = f{x)'^x^~^ dx. Once a formal expansion of the energy is achieved, one can easily demonstrate that a ^ 

One can then show that there is an upper bound on the I^-dimensional energy eq. ([5]) of the form 

E < E{a) = eo + ei + 62 + 63 (8) 



where the constants Cg are given by 



^2 roo 



(9) 
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It is clear that eo, 62, and 63 are all positive numbers, but ei = ki (70 — A^) may become negative if A^ > 70- For example, 

2D/2 



for a Gaussian wavefunction f{x) = e ^\Md = Vd/1+\ ^ ^nd one can deduce the following constants: 



Dn 
27^' 



27r 2 



^2 



I^n^ A^ 

4 7r 2 '?/;o 

A^ 

3^/4 • 



(10) 



It should be stressed that the variational approach only provides a ceiling on the ground state of the energy. If this upper bound 
— the lowest variational energy — is negative for a 7^ (since, for the nearly uniform state, the energy is positive and tends to 
zero as a ^ 0), then the uniform state is known to be unstable and it is most probable that the impurity field has self localized. 
However, should the upper energy bound be non-negative, one cannot determine from this whether or not self localization is 
likely to occur. 

Regardless, the energy expression eq. dU) in terms of the variational parameter a does provide some guidance as to when a 
localized impurity solution is likely to exist in 1, 2, and 3 spatial dimensions. The key here is the observation that when ei > 0, 
the function ([5]) increases monotonically, implying that the energy is minimized when a ^ 0. We will now examine the situation 
of one, two, and three dimensional space in turn. 



A. Self localization in ID 



Applying this in one dimension is straightforward: since when D = 1 the dominant term of E{a), i.e. eq. ([8]), at small a is 
ei a, when ei is negative, i.e. A^ > 70, a supercritical transition occurs from a homogeneous state to a localized impurity state. 



B. Self localization in 2D 

For D = 2, the first and second terms of E{a) in eq. ([8]) are of the same order. As occurs fovD = 1, there is a second- 
order transition towards a localized impurity state if A^ > 70 + , with ki = AC2 A/2 hoa^r^^\ . The instability of the 

homogeneous state is shifted from the bulk condition (A > 70). This shift has a simple interpretation: the presence of a 1/m 
factor means that this term comes from the kinetic energy of the impurity / 2^ | VxP, so the shift is created by the curvature a 

f°° f' ix)'^ X dx 

of the localized structure — it is a kind of surface tension. Furthermore, ?oo p/ — 1— bears the hallmark of a surface tension 

Jo fixrxdx 

effect: that is the ratio of interface energy (the gradient term) to a bulk energy. 

In the case of the Gaussian trial function one finds that a localized solution exists if A > Ac = \JlQ^r^^- These criteria are 

close to those numerically observed in two-dimensional space. Indeed, we have simulated numerically eqs © and ^ for the 
case of a single impurity {N = 1) with 70 = 0, m = 1, and = 40.96 in a 64^ box in which the total number of particles in 
the condensate was J Itpl'^cPx = 4096. Below we present our findings for two different sets of initial states: i) uniform miscible 
state, i.e. ^ 1 and x ^ 0-1. both real (here ^ C means C plus a small fluctuation or noise); and ii) = I, and x is set as the 
Gaussian Ansatz ^ where the value of a is that which minimizes the energy E{a) in 2D. By the variational criteria, a localized 

solution should exist if A > . ^ 0.3917. 

For the initial condition i), we find numerically that if A < 0.5 no instability will develop because of the finite size of the 
system. Indeed, if A < ^ 0.491, the system is not unstable for the longest wavelength allowed by the box (a wavenumber 

^). However, as soon as we get a localized solution, e.g. A = 0.5, it is possible to decrease A to below Ac = 0.39 and retain 
the localization. Indeed, we can bring the wavelength down to as low as A 0.35 and still keep the nonlinear solution. Below 
A 0.35 it is not possible to distinguish a real localized structure from oscillatory structures. 

For A = 0.35 to 0.55 an oscillatory pulse is observed — a kind of breathing. This is perhaps due to an excess of energy, 
and the localized solutions must radiate to evacuate it. The nature of the oscillation should follow from a similar Ansatz to 
eqs (I6I7I) but allow the phase of the impurity field to be a space- and time-dependent function. For instance, taking x(^, t) = 
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V Cd^^d r)e*^^*^^^ and introducing this into the Lagrangian of the system (which is — / | {xdtX~xdtX)d^ ^ minus the 
energy eq. ([5])), for the Gaussian trial function one arrives at the following Lagrangian: p, p) = —E{a) — ^2 + f ) • 

The Euler-Lagrange conditions lead to the "momentum" equation p = — Finally, this yields the dynamics that the "breath- 
ing" must satisfy: 

VA — Ac and the breathing frequency is oobreathing ^ 

♦ 

FIG. 1: (Color online) Numerical simulation (see Appendix B for details) of eqs ^ and 0]) for a single impurity (A^ = 1) with 70 = in a 
64 X 64units^ periodic plane. Here J {tpl'^d^x — 4096 and J \xkfd^x — 40.96. The plots represent |xi|^ and the colormap is the same for 
all four images. In aj, A = 0.385 (no localized solution is observed, but because of the recurrence phenomenon discussed in the text, after a 
time a kind of pulse is observed, which then dissipates, only to eventually form again); in b), A = 0.39; in c), A = 0.395; and in d), A = 0.4. 
All these solutions are dynamically obtained from the same initial state i) and, as described in the main text. We start with A = 0.5 and, after 
the spontaneous formation of a localized solution, we decrease A slowly, i.e. A varies in time every 1000 time units by A ^ A — AA where 
AA = 0.005. 

For the initial condition ii), valid only if A > Ac, one notices that the Gaussian trial function is very close to the exact stationary 
solution, and the pulsations and breathing described previously are less marked. From an initial condition of A > Ac, we can 
decrease A to below Ac, retaining the localization. The oscillations, however, become very important and, as earlier, it seems 
that rather than a stationary state, a more complex oscillatory behavior is displayed. This dynamic solution seems to be very 
robust as we can observe it even for very low A. The pulse forms, disperses through the system, and then after some (possibly 
long) time the pulse forms again, in a kind of recursive, and not necessarily periodic dynamic. However, this dynamic behavior 
dominated by pulses seems to be quite different from the steady ground state sketched by our variational theory. 

Nevertheless, we notice that the numerically observed range over which steady-state localized solutions exist is very close to 
that expected via the variational analysis. The difference probably arises because we are simulating the Hamiltonian evolutions 
eqs ^ and (O, and not simply finding the minimum energy of the Hamiltonian eq. ([T]). 

C. Self localization in 3D 

For D = 3, finding the self-localization conditions is more subtle than in lower dimensions. As discussed above, the vari- 
ational argument shows that a localized solution exists if the energy E in eq. ([8]) is negative for some critical value, ac. A 
necessary (but not sufficient) condition for a negative-energy ground state in 3 D is ei < 0, i.e. A^ > 70. 

We can make a more precise estimate as follows: Because the energy turns negative for a > 0, and because in general the 
energy expansion grows as near a = (since eo > in eq. ([5])), it is sufficient (though not necessary) that E{a) have a 
minimum at a = o^c and that this minimum reach the horizontal axis E{ac) = 0. In equations, this reads as E{ac) = and 




6 



E'{ac) = 0, which describe a critical Hne in spatial parameters, separating the region where a localized structure (with great 
probability) exists from that where there is no such certainty. This line may be written parametrically as 

3/4 / — 

6/3^ + 2/33 = !2l2_ & 7/36 + 3/32 + !iV£3 ^ ^^^^ 

where f3 = ac (63/62)^^^ is the parametrization. The parametric curve may be approximated in the large and small f3 limit 
respectively by: 

7 6/7 1/7 

= - '"ge/? (12) 

o 2/3 1/3 

^1 = - %2/l' if/3«l- (13) 

Let us look briefly at the case 70 = 0. Here, one may express eq. (fTTt as a closed parametric expression for the variables 
(m A , ). The values of the e's in eq. (fTOl) for the Gaussian trial function yield a bound that is about 35% higher than the 
numerical result flTI] , as plotted in Fig. [2 Moreover, for a Gaussian trial function the asymptotic formulas ([T2b and ilSh take on 
the simplified expressions 

9^\^/^ 1 1 



mA = — . = 7.2907 ^ = for (n^V^oM <Cl, (14) 

77/10^3/5 ;l 1 
= 017/2Q 7 T~ryni ^''^^^^l / / A2/5 for(n*V^o/m) > 1. (15) 




FIG. 2: (Color online) Parametric plot of the critical transition line for 70 = in three spatial dimensions. The curve /) is the parametric curve 
fTTI) : the curve //) is the asymptotic behavior eq. ([T4l) for n^^g/m ^ 1; and is the asymptotic behavior eq. (fTSl ) forn^^o/m ^ 1. The 
curve iv) represents the Cucchietti-Timmermans condition fllll : -^X^n^il^Qm > 4.7, that is Act = ^nj^^gm ' Fi^^HY' the inset represents 
the same curves, but in log-log to cover a larger range of values. 



V. MEDIATED ATTRACTION OF IMPURITY FIELDS DUE TO THE CONDENSATE 



In this section, we look at how N localized dilute impurity fields interact and determine the role of the modified background 
condensate in this interaction. We will demonstrate that the interaction between impurity fields has an attractive tail, and that 
this is a result of mediation by the condensate, combined with the hard-core repulsion arising from the positive scattering length 
between impurity fields. 

Consider N impurity fields that are self localized (having satisfied the conditions given in the previous section) and that 
weakly modify the initially uniform condensate, i.e. 



(16) 
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with ipo assumed constant. This assumption holds if 1) there is only weak localization of the impurity or 2) the impurity is 
distant enough. The second case is particularly relevant here as we intend to derive the nature of the interaction at long range. In 
this approximation, using eq. © the condensate wave function is given by (here | Xfc P is assumed to be of order e) 



1 ^ 



(17) 



k=l 



This linear equation (known as the screened Poisson equation) may be solved with the aid of a Green's function in I^-spatial 
dimensions, 

-AG^^\x - x') + {2^o)^G^^\x - x') = CdS^^\x - x'), 
where Cd the surface of a unit sphere. The solution can be written explicitly as 



G^''\x-x') = 



( -^e-2^ol^-^'l inD = 1 

Ko{2ijo\x-x'\) inD = 2 

-2iPq\x-x'\ . 

I .1 — inD = 3 

\x—x' I 



(18) 



where Kq is the modified Bessel function. 

Solving eq. dTTl) yields the following modification (which is of order e) to the condensate field: 



2Xipo 
' Cd 



(G(^\x-x')Y,\Xk{x')?d^x'. 



(19) 



Under the assumption of a weakly modified condensate, the total energy (on the order of e^) is given by 



E 



. / N N N N 

y /c=l /c=l k=l j<k 



dx. (20) 



The consistent scaling is achieved by assuming |Vx/cP ^ (which plays a role in the self energy as discussed below) and 

\Xk? ^e. 

One can multiply eq. ([TtI) by V^i and integrate over the total volume to obtain 

j (^^IV^Z-il" + 2V^2|^i|2 + A^/'o^/'i ^ \Xk\^ dx = 0. 
Using this identity we can eliminate the terms in eq. ([2Ql) that do not depend explicitly on the impurity fields to arrive at 

. I N AT AT AT \ 

E= / A^o^i^|xfeP + ^— |Vxfc|2 + |5:iXfc|" + 7Elx^-|'IXfcl'H''- 

\ fe=l /e=l j<k ) 



(21) 



Utilizing the solution for V^i(r), i.e. eq. ([T9l) , the energy of the system of N impurity fields under the assumption of a weakly 
modified condensate can be written as 

= E/ (1^5^''\x-x')-'^-^G^^\x-x')^ \xM')?\Xk{x)\U^xd^x' + E, ill) 
where the self-interaction energy of the impurity fields is given by 

= (i'^^^'' + T'^^0 ^"'"^ lG'-^Hx-x')\xk{x')\'\xk{x)\'dxdx' 

(This self-interaction term Eq was previously derived in ifTTI] .) 

Let us consider localized impurity fields, a\xk — Xi\ ^ 1. One may approximate the k-th impurity field by 

\Xk{x)\^ = nkS^^\x-Xk) 



(23) 
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where Xk is the position of the k-th impurity. In this case E = ^ ^i^k ~ ^k\)^i'^k where the interaction potential 

between the two impurity fields is given by 

U(\x, -Xk\)= i8^''\x, - Xk) - —^G^'^Hx, - Xk). (24) 

Note that the self-interaction energy Eq naturally diverges because of the artificial singular structure \xk{x)\'^ = nk5^^\x — Xk) 
assumed in the impurity field. 

In eq. (l24l) . the first term on the RHS is, on its own, merely a crude estimate. The full energy expression eq. (l24l) points to the 
existence of a bound state but with an equilibrium distance of zero. In realistic systems, the equilibrium distance is dependent 
on the impurity's localization, which is on the order of the size of the structure. 

For example, consider the trial function eq. © of the form Xk(x) = V^yc^^^^""^'^"^'"'^- Introducing this into the 
interaction energy eq. (l22l) , one arrives (in 3D) at 

U{\x, -Xk\)= 7^e-^l---^l^ - ^ , . (25) 



Xk\ 



The second term is estimated in the limit a ^ 2?/^o, i-e. where the structure is more localized than the Yukawa interaction range 
of l/(2?/^o)- The equilibrium distance comes from this energy. 



VI. N IMPURITY FIELDS - CRYSTALLIZATION 



Here we will discuss a system composed of N interacting impurity fields inside a larger condensate. In section |Vl we showed 
that over a certain range of parameters the tunable interaction between localized impurity fields has a hard core and an attractive 
tail, which is a type of interaction that has been studied in the context of many diverse physical systems (for Yukawa-type 
attractive interactions see e.g. ifTsIl ). We show below that, depending on the interaction between condensate and impurity, there 
are two possible regimes in which such impurity fields crystallize: the first where the impurity fields remain immersed in the 
condensate, and the second where they phase separate in the form of a crystal within a bubble within the condensate. 

The phase diagram in Fig. [3] shows four distinct regimes. In phase /, the impurity fields are miscible with each other and with 
the condensate. In phase //, the impurity fields are miscible with each other but phase separate as a bubble from the condensate. 
In phase ///, localized impurity fields stay miscible within the condensate and can crystallize if system is not dilute; in phase IV 
the impurity fields phase separate from the condensate and form a crystalline structure within a bubble 

. / (]\f — \\ -y-|--yQ I I 

In an infinite system, the critical lines for the equilibrium phase are 7 = 70 and A = — as discussed in section Hill 

In the following subsections we shall discuss the different phases. Although in a finite system these transition lines are shifted 
because the lowest mode (where the instability typically arises) is nonzero, for simplicity, the transition lines discussed will be 
those for the case of an infinite system. However, it should be kept in mind that corrections maybe computed (See Appendix A 
for discussion). 



A. Phases I Sell 



In the case 7 < 70 and A < ^ , the condensate-impurities system finds its minimum energy in a homogeneous 

state, which is the completely miscible state. This is phase / of Fig. [3l In phase //, 7 < 70 and A > and 
the condensate phase separates from the impurity fields, which remain together and miscible with each other in a bubble. No 
crystallization occurs in either phase. 



B. Phase /// 



In this case, 7 > 70 and A < JlS^^'^^ . In phase ///, the interaction energy eq. (1221) causes the impurity fields to attract. 
However, because a repulsive hard core exists at close range, one expects there to be an equilibrium distance. Although we 
expect the system to crystallize, the complete picture is not so simple because there is also a density parameter that plays a 
fundamental role. Indeed, as already seen in section |lVla~^ is roughly the size of a localized structure (see section Hvl). so then 
the "diluteness" of the impurities in the condensate is measured by ^^^p"' ^^^^ parameter is close to unity, as is true in typical 
solids, one gets a crystalline phase (see Fig. (4]). However, if this phase's "diluteness" becomes small, one would expect to see 
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FIG. 3: (Color online) Phase diagram for = 6 impurity fields with 70 = 1. The vertical line i) represents the border line 7 = 70 = 1 while 

, / ( iV — 1) 7+70 

the curve ii) represents the critical line — . The accompanying pairs of pictures (of which the left represents the condensate density, 

and the right plots J2k=i \Xkf) ^^e numerical simulations in a 64 x 64 unit periodic plane with uniform (plus a small fluctuation) 
initial conditions for the condensate and the impurity fields such that f \ip\'^(Px = 4096 and f \xk\'^d'^x — 40.96. Explicitly, we have 
A = 0.25 and 7 = 0.5 in /; A = 2 and 7 = 0.5 in //; A = 0.25 and 7 = 2 in ///; and A = 2 and 7 = 2 in IV. See Appendix B for details 
of numerical tools. 



liquid- or gas-like behavior, as occurs in molecular dynamics simulations [18]. These states are roughly observed in numerical 
simulations (see Fig. O. 

C. Phase IV 



In this phase, 7 > 70 and A > — . The condensate phase separates from the impurity fields, and the impurity 

fields also phase separate from each other (unlike in Phase II). This results in the impurity fields gathering into a "bubble" inside 
of which the condensate density is essentially zero. Within this bubble, the absence of mediating condensate means that the 
impurity-impurity interaction is no longer Yukawa-type (i.e. ?/^o ^ in eq. (fTSl)), a special case which will be discussed in the 
next section. The condensate basically acts as a container by causing a stronger attraction between any impurity field that strays 
from the bubble and the remaining ensemble of impurity fields. This inward "pressure" from the condensate and the close-range 
hard-core repulsion bring about crystallization (see Fig. [6|. 

VII. ANALYTICAL TREATMENT OF PHASE IVi THE CASE A = 

In phase IV described above, the system is at a special limit when the condensate is not coupled to the impurity fields. In 
such a case, the induced attractive interaction described earlier is no longer valid. Therefore, the case A = deserves special 
attention. The condensate can be described by the pure nonlinear Schrodinger equation, which is quite well understood. The 
impurity fields, however, evolve and we shall now consider this evolution. Impurity fields that overlap store potential energy and 
will tend to repel each other; more so as 7 increases. If the impurity fields \xk\^ are localized about Rk with a size 3 then, in the 
limit 7 ^ 7o, the potential energy is dominated by 

|E/ \Ux)\'\Xk{x)\'d^x. 

Minimizing the overlap will minimize the potential energy. In some sense, the energy may be approximated by the superposition 
of two-body interactions £ = ^ Yl^i^k ~ with U (d) repulsive and Gaussian (see below). However, in a finite box 

there is a limit to how far away from each other the impurities can move. A kind of close-packing argument may be invoked to 
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FIG. 4: (Color online) Numerical simulation of eqs ^ and ^ for N = 6 impurity fields in a 64 x 64 units^ periodic plane. Here the initial 
state is given by uniform (plus a small fluctuation) wave functions for the condensate and the impurity fields such that f \ip\'^(fx = 4096 and 
/ \Xk\'^d'^^ — 40.96, moreover m = 1, A = 0.2, 7 = 1 and 70 = 0. Figure a) plots the condensate density |^|^, and b) plots J2k^i IXfc|^- 
Figures /) through vi) each plot a different impurity field |xfe See supplementary material for a movie of crystallization Ii l9i1 . 



find the crystal structure that minimizes the energy of the impurity ensemble, sustained only by the external pressure from the 
condensate at the boundaries. 

The periodic case, i.e. when all TV impurity fields have the same number of particles n*, provides an interesting, solvable 
example. One can use the minimization approach to determine the size ^ as a function of the large parameter 7. Only the 
nearest neighbor(s) affect the interaction energy. To begin, let us assume that the impurity field is a compact function (i.e., it 
vanishes exactly outside the ball of radius 6) and that there is no overlap. In this case, one can calculate the minimum of the 
total energy. We shall assume 70 = in the following and note that there is no overlap, which means no interaction energy 
2 Xli//c / \Xi{^)\'^\Xk{x)\'^d^x- The total energy can be written as 



H 



N ^ 



\Vxk?dx. 



The Euler-Lagrange conditions for an extreme of this energy leads to the Hemholtz equation 

- -^^Xk = exk 
2m 



(26) 



(27) 

inside a ball Vs{Rk) and with a Dirichlet boundary condition on dVs{Rk) II2QIi . 

From here, one can proceed using perturbation to determine the effect of a slight overlap. The variational parameter S will 
determine the optimal configuration that balances the kinetic energy of a pulse and the interaction energy due to a small overlap. 

As an example, let us consider the one-dimensional case in a periodic domain of length L. Let Xs{x) = ^/n^f{x — sd) with 
d = L/A/', the ground state. The energy is dependent on the number of impurity fields (which is the same as the number of sites) 
such that 



E[f] = 



The energy may be estimated explicitly by using 



cos(7r^/(2(5)), 



dx 



(28) 



(29) 
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FIG. 5: (Color online) Numerical simulation of eqs Q and 0]) for N = 6 impurity fields in a 64 x 64 unit periodic plane. Here, as in 
previous figure, the initial state is given by uniform (plus a small fluctuation) wave functions for the condensate and the impurity fields such 
that / l^l^fd^x = 4096 and f \xkfd^x = 40.96, and m = 1, 7 = 1 and 70 = 0. A increases from left to right. The upper row plots the 
condensate density \ip\'^ while the lower row plots J2k^i IXfc The corresponding values of A are as follows: A = 0.4 in i) and ii) ; A = 0.5 in 
Hi) and /vj; A = 0.7 in vj and v/j; and A = 1 in vii) and viii). The size of the impurity fields decreases as A increases due to the self-interaction 
term seen in eq. ( [23] ). (In the latter case, notice that only five of the impurity fields interact strongly and ultimately form a bound state; the 
remainder are far enough so that the exponentially weak interaction is negligible.) The case A = 0.2 is plotted in a) and b) of Fig. |4] See 
supplementary material for a movie of crystallization fioll . 




FIG. 6: (Color online) Numerical simulation of eqs ^ and ([S]) for = 36 impurity fields in a 64 x 64 units^ periodic plane. The initial 
state is given by uniform (plus a small fluctuation) wave functions for the condensate and the impurity fields such that J = 4096 and 

/ \xk\'^d'^x = 40.96, A = 2, 7 = 1.5 and 70 = 0. Plot a) depicts the condensate density |^|^ and plot b) depicts J2k^i IXfcl^. 



the solution of the ID Hemholtz equation (|27l) . This yields 



E = iVn 4 [\-\rs{x)\^^jnJs{x)\f5{x-df^fs{x-df)]dx 
2 J-s m 



n*7 



(3 S sm{7rd/S) - tt (d - 2 J) (2 + cos{7rd/S))) 



The dimensionless parameter ( = ^ that minimizes this energy satisfies the implicit equation 



mn^^d 



(-4 ^ + 4 C - (2 ^ + cos(C) + (-3 + 2 ^ C - C') sin(C)) 
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d 



-< >- 




FIG. 7: (Color online) Sketch of the minimum energy (overlap) configuration for two impurity fields. The size of the localized fields is 5 and 
they are separated by d. 



which in the large 7 limit (C ^ 27r) gives 
Thus, as expected, (5 ^ (i/2 as 7 ^ 00. In a similar approximation, the energy becomes 



d 
2 



{mn^jd) 4 + . . . 



E = Nn^ 



2md^ 



1/4 



(30) 



The compact solution f5{z) is only an approximation since the real ground state does not vanish entirely at x = ±6 even 
though the impurity field is exponentially small. The nonlinear structure of the boundary layer may be approximated by the 
ordinary differential equation that satisfies the function / : 



1 

2m 



f" + inkf{x)U{x - df + f{x + df) = ef{x). 



We are interested in the region where x ^ 5, so one may neglect the f{x + d^, which is an exponentially small term. We 
approximate the second term by fs{x — d)'^ from eq. ([29l) , which is a given and known function. 

The remaining linear ordinary differential equation may be solved in the WKB (large 7) limit, f{x) ^ e~^^^^\ so that 

^S\xf ^ rikfsix - d)\ Therefore, 

b{x) ^ -\/2mnkOsm{ — ). 

TT ZO 



Finally, near x ^ S ^ d/2 (sls shown previously S ^ d/2), we find 



The domination of the interaction U (d) by this Gaussian tail shows that impurity fields can be free of each other's influence 
with sufficient separation. Also, when the separation is sufficiently small, the hard-core repulsion term comes into play, making 
the impurities impenetrable (see Fig. [8]). 

In conclusion, a crystal phase is possible whenever the condensate presence is negligible and the impurity fields are immisci- 
ble. The crystal is composed of distinct localized impurity fields which are repulsive and impenetrable. Therefore, the crystal 
lattice is a result of the effect of the boundaries: if V is the volume then the mean separation distance would be ~ {V/Ny^^ 
which determines the localized structure size. This crystal will not exist without the presence of boundaries or a containing pres- 
sure. The condensate acts naturally as this pressure as seen in Fig. [6l By construction, it seems plausible that these crystallites 
(like those of Figl3]-/y, Fig. ^viii), or Fig. [6]) possess a surface energy. In general, the energy of the system in phase IV is made 
up of a bulk contribution from the condensate component (i/j ^ i/jq and Xk ^ 0) set to zero in the energy given in eq. (0)), 
the energy of the crystal phase (ip ^ and Xk ^ crystal structure built in the way explained in this section), and an interphase 
energy, which we have not computed. However, by general arguments one can expect that the ratio of this interphase energy to 
the bulk energy from the crystal gives a critical radius of crystallite nucleation. It also says that the transition between phases /// 
and IV is of first order. However, we were unable to observe a hysteresis effect (in Fig. O when the crystallite is formed and one 
decreases A through the boundary. This implies the interface energy is weak. 
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FIG. 8: (Color online) Numerical simulation of eq. 0]) for = 36 impurity fields with A = in a 128 x 128 units^ periodic plane. Here 
the initial condition is an uniform miscible state plus small fluctuations, the number of particles are given by f \xk\'^d'^x = 200, 7 = 1, and 
70 = 0. The plots represent the instant t = 1700 units : a) shows J2k^i IXfc 1^^ i) H) plot two of the 36 distinct impurity fields. Note 
that the system evolves through a Hamiltonian (conservative) dynamics though the relaxation to an "equilibrium" takes a long time. Indeed, at 
the moment depicted in the plots above, the system has yet to reach such an equilibrium, so each impurity field contains many distinct localized 
structures. We remind the reader that while the analysis done above is in ID, the numerical simulations are in 2D. 




FIG. 9: (Color online) Numerical simulation of eqns. JJ]) and 0]) for {N = 6) impurity fields with 70 = 1 in a 64 x 64 units^ periodic 
plane where / |^|^(i^£c = 4096 and / \xk\'^d^x — 40.96. We begin with a crystallite initial state that comes from the same simulation as in 

Fig. [3l-phase IV, and subsequently we slowly decrease A from A = 2 to below the critical value — ^JllJQ ^ 1.354. The plots 

represent \Xk\^ and the colormap is the same for all four images. The parameters are A = 1.5 at the stage shown in a), A = 1.3 in b), 
\ — 1.25 in cj, and then we increase back to A = 1.3 in d). We do not observe a hysteresis effect. 



VIII. NONCLASSICAL ROTATIONAL INERTIA 



In this section, we will investigate the crystaline phases in phases /// and IV retain their superfluid behavior, thus exhibiting 
properties reminiscent of a supersolid dill 0, |23|] . Inspired by Leggett's seminal work ll24l] . consider a condensate-impurity 
system confined within an annulus of total length L and total volume V = SL. Consider this system to be rotating uniformly 
about the annulus's primary axis of rotation, ei (the direction associated with the coordinate xi which spans [0,L]), a rotation 
induced by the boundary condition = L) = t/j^xi = 0)e*"° where ao = is dimensionless, and X/c(^i = L) = 

Xk{xi = 0)e*"^"^p where aimp = . Then the system will possess an energy (for ao <C 1 and aimp <C 1) 
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where Eq is the ground state energy and AS = ^I^^^uo'^, being the effective and observable (measurable) moment of inertia 
tensor around the ei axis. 

The relative deviations of this tensor with respect to the rigid body rotation, I^^ = pL'^V where p is the total mass density, 
is called the nonclassical rotational inertia (NCRI) fraction of the annulus. NCRI is seen as a signature of a superfluid response 
of a system. 

Let i/j^ and Xk the ground state real wave functions of the Bose field and the impurity fields respectively. Under rotation, 
the phases of these fields are no longer uniform, so the rotating wave functions are of the form ilj^{x)e'^^ and Xki^)^^^'' (the 
only change comes from the phase to the lowest order) and the increase of energy is 



A£ = ^ I [ po{x) |V(/)|^ + V -pk{x)\\/M^ 1 dx 




m 



where po{x) = |?/^^(x)p and pk{x) = \x^{x)\'^ are the non-uniform density of the respective ground states, and the phases (j) 
and (j)k satisfy the boundary conditions explained below. That A£ is a quadratic form in cj, i.e. in a, is quite evident since (j) and 
are proportional to ao and aimp because of the boundary conditions (see below). The explicit form of I^^ ^ and the prefactor, 
on the other hand, need longer consideration. 

The minimization of AS leads to the equations and boundary conditions (recalling the coordinates other than xi are periodic) 

V-(po(ic)V0(x)) =0 eL^ & 0(xi = L) = ^(xi = 0) + ao (31) 
V • {pk{x) VMx)) =0 eL^ & Mxi =L) = (})k{xi = 0) + aimp- (32) 

Eqs (|3T1) and (l32l) may be solved using the method called homogenization [25]. This method splits cleanly the large (system 
size L) and small (impurity size) scales and provides effective average quantities. Consider for example eq. (|3T1) , taking 
(j){x) = + i>{x) where ^{x) is periodic in . The periodic function ^{x) = aoK{x) satisfies 

V . {po{x) VK{x)) + ei • Vpq{x) = 

and, after substituting this back into Af , one finally obtains 



where 



^0=^1^ Po{x){l - {yK f)d^x (33) 

and has units of a mass density. In fact, it is the superfluid density. Notice that by replacing one can identify the effective 
moment of inertia I^^^ = q^^L'^V, which is the moment of inertia of an annulus with effective density mass g^^ . It is not 
possible to obtain a closed expression for g^^ in terms of the local density Pq{x) (K implicity depends on po). However, it is 
possible to show that I^^^ < I^^ , and hence the system displays a NCRI. In one-dimensional space, eq. ([3T1) may be solved 
exactly, with a direct calculation leading to the formula 

\V Jv Po{x) J 

due to Leggett iH. 

Finally, to take into account the contributions from the condensate and impurity field, one need simply add the corresponding 
superfluid densities; and because the impurity fields are localized, their wave functions decay fast in space so the contribution to 
the superfluid density from the impurity fields becomes negligible. Therefore, perhaps unsurprisingly, the superfluid density is 
dominated by the condensate density. 



IX. DISCUSSION 



In this paper, we analyzed how the condensate-mediated interactions between self localized impurity fields can lead to pattern 
formation. We restricted ourselves to the case where all impurity-impurity interactions are governed by the same coupling 
constant. In ultracold trapped atomic gas experiments, impurity fields may be generated by populating different atomic levels 
and/or by making use of different atomic species/isotopes (for review of current experimental progress on multicomponent 
condensates see section IX of IH] ). In these cases, unlike pairs of impurity field would have different coupling constants. In this 
situation, we would expect there to be phases similar to those described in sectionjVll but that the broken translational symmetry 
of the Hamiltonian will lead to nonperiodic structures, rather than periodic order. 
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APPENDIX A: DYNAMIC INSTABILITIES OF THE MISCIBLE STATE: THE BOGOLIUBOV SPECTRA AND FINITE SIZE 

EFFECTS 

In this Appendix we work out the dynamic instabihty conditions of the uniform miscible state, and show that these conditions 
are important when we take into account the fact that the system is of finite size. Consider the case of uniform ground states, 



iQnt 



where l^o = IV'oP + AEf=i \x°\\ and Xj = x^e"^"^* where % = A|VoP + 7o|x°P +lEk^j 

1 BogoHubov dispersion relations [13] , 



+ 1 modes with wavenumber k around this ground state would obey the N 



The 



where s runs from 1 to AT + 1 and e^*^ is the s^^ eigenvalue of the matrix 



/I 

1/m 



Mb = 



1/m 



AlxiPM 7o|x?lV^ tIXiPM 
MX2?/m 7\X2?/m 7o\X2?/m 



AlV^ol 



l\Xi?/m 



7|XilV^ 
l\X2?/m 



: 1/m A|x^_ilV^ : 7o|x^-ilV^ tIx^-iIV^ 

VO l/m/ \ X\x%\Vm j\x%\Vm tIx'^IV^ 7o|x'^lV^ 

(Al) 

From the construction of the matrix it is easy to show that if M, eq. dH), is positive semidefinite then Mb^^ also positive 
semidefinite. To do so, we note first that the sum of two positive matrices is a positive matrix and that the first matrix on the 
RHS of eq. (lAll) is naturally positive; then, for the second matrix on the RHS of eq. (lAll) . we recall that a matrix is positive 
semidefinite if the determinants of all upper left submatrices are non-negative and note that the determinants of all the upper 
left submatrices of Mb differ from the corresponding submatrices of M only by a factor of |?/;oP or |Xj P, both of which are 
positive. 

The eigenvalues of this second matrix play an important role in the dynamics because they are they charaterize the instabilities 
at long wavelengths. In the convenient case where all |Xj P = Ix^P' the eigenvalues are (70 — 7)|x^P/^' which is A" — 1 
degenerate, and 



i (iV^oP + {{N - 1)7 + 7o)|x'lV^ ± V(-|V^oP + {{N - 1)7 + 7o)|x^l V^)' + mX^l^olW/rn 



(A2) 



V^(Ar-l)7+7o 



/N 



, which are the same as the energetic 



We can see that a modulational instability occurs if either 7 > 70 or A > 
instability conditions described in section HIH 

Therefore, the dynamic instability conditions (from the Bogoliubov analysis) and the energetic instability conditions are 
identical in the infinite system. Note that the values of the densities of the condensate \ipo\'^ and impurity fields |x^P are absent 
in these conditions since they only play a role in the time and length scales involved. 

However, in a finite system — take for instance a periodic box of size L these transition lines are shifted because the lowest 
mode (where the instability typically arises) is nonzero. Indeed, the growth instability rate becomes positive if one of the eigen- 
frequencies uol < 0. That is, the critical line is shifted by j^. So, assuming |Xj P = |x^P as before, the degenerate eigenvalue 
(70 — 7)|x^P (responsible for the impurity immiscibility) becomes negative when 70 + l^Jx^i"^ ^ 7 (instead of 70 < 7). 

One can also show the instability condition responsible for the condensate-impurity immiscibility is shifted as can been seen 
in the large L limit: 



A> 



V(A^-l)7 + 7o , |V^oP + ((A^-l)7 + 7o)|x'P 



N 



• 0{1/L^) 



(A3) 



where the latter approximation is valid in the limit |x | ^ IV^o| • One notices that this second condition is shifted by a 
smaller amount as the number of impurities A^ becomes larger because of the ^/W(W^^T) in the denominator. 



APPENDIX B: NUMERICAL TOOLS 



We used a Gauss-Seidel Crank-Nicholson finite difference method to integrate the Hamiltonian coupled equations (I2I3I) . The 



scheme is as follows: let U{t) = Xi(^) • • • Xa^(^)} be the fields at time t, and let us write eqs (I2I3I) as 



du _ 

dt 



F[Uit)] 
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where F[U] is a. nonlinear Hermitian operator whose definition is clear by comparison with eqs i2\3i . Thus, one time step from 
t to t -\- dt is 

U(t + dt) - ^F[U(t + dt)] = U{t) + ^F[U{t)]. 

For small dt, this (classical) numerical schema yields ip{t -\- dt), exact to at least the second order in dt and has the advantage 
that, formally, the norm and energy are exactly conserved. However, the price is that U{t -\- dt) involves calculation from a 
complex nonlinear equation. Next, to obtain U{t-\- dt), one iterates the mapping Un-\-i = {U (t) + ^^[^ (^)]) + ^F[un] where 
{U (t) + ^F[/7 {t)]) is constant, Un^o = U{t), and Un^oo = U{t -\- dt). Convergence is expected after a moderate number of 
iterations, 6 in practice, using a time step dt = 0.01 units. Under this condition, the norm and the energy of the solution are 
conserved in time, deviating by less than one part per 10^ per time step and less than one part per 10^ per unit time. 

For most of the figures, the numerical simulation of eq. ^ was performed in a 64 x 64 units^ periodic plane for = 6 
impurity fields. Unless otherwise stated, the simulations typically use / \ilj\'^d'^x = 4096 and / \xk\'^d'^x = 40.96. 
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